function [ksi]=bound_coeffs(mesh,sp)

%Calculates boundary condition coefficients based on reflectivity.
%Performs numerical integration of Fresnel law.

nm=mesh.ri;
n0=ones(length(nm),1);
Ca=zeros(length(nm),1);
Ca(:)=cos(asin(n0./nm));       %w relating to Critical Angle


%only calculates integral once for each unique RI.
ri=unique(mesh.ri);

Rn1=num_int(ri,n0,Ca,ones(size(ri)));
Rn1=reshape(Rn1,length(ri),15);

Rn=zeros(length(nm),15);
for i=1:length(ri)
    ind=find(nm==ri(i));
    for j=1:15
	Rn(ind,j)=Rn1(i,j);
    end
end

Rn=reshape(Rn,length(nm),15);



%Create output matrices for boundary coefficent values

if sp==1
    
    BC.A=-Rn(:,2);
    BC.B=3*Rn(:,3);
    
    ksi=[BC.A BC.B];

    clear BC

elseif sp==3
    
    % A1,A2
    BC.A(:,1)=-Rn(:,2);
    BC.A(:,2)=(-9/4)*Rn(:,2)+(15/2)*Rn(:,4)-(25/4)*Rn(:,6);
    
    %B1,B2
    BC.B(:,1)=3*Rn(:,3);
    BC.B(:,2)=(63/4)*Rn(:,3)-(105/2)*Rn(:,5)+(175/4)*Rn(:,7);
    
    %C1,C2
    BC.C(:,1)=(-3/2)*Rn(:,2)+(5/2)*Rn(:,4);
    BC.C(:,2)=(-3/2)*Rn(:,2)+(5/2)*Rn(:,4);
    
    %D1,D2
    BC.D(:,1)=(3/2)*Rn(:,3)-(5/2)*Rn(:,5);
    BC.D(:,2)=(3/2)*Rn(:,3)-(5/2)*Rn(:,5);
    
    ksi = [BC.A BC.B BC.C BC.D];

    clear BC

elseif sp==5
    %A1, A2, A3
    BC.A(:,1)=-Rn(:,2);
    BC.A(:,2)=(-9/4)*Rn(:,2)+(15/2)*Rn(:,4)-(25/4)*Rn(:,6);
    BC.A(:,3)=(-225/64)*Rn(:,2)+(525/16)*Rn(:,4)-(3395/32)*Rn(:,6)+...
        (2205/16)*Rn(:,8)-(3969/64)*Rn(:,10);
    
    %B1, B2, B3
    BC.B(:,1)=3*Rn(:,3);
    BC.B(:,2)=(63/4)*Rn(:,3)-(105/2)*Rn(:,5)+(175/4)*Rn(:,7);
    BC.B(:,3)=(2475/64)*Rn(:,3)-(5775/16)*Rn(:,5)+(37345/32)*Rn(:,7)-...
        (24255/16)*Rn(:,9)+(43659/64)*Rn(:,11);
    
    %C1, C2,C3
    BC.C(:,1)=(-3/2)*Rn(:,2)+(5/2)*Rn(:,4);
    BC.C(:,2)=(-3/2)*Rn(:,2)+(5/2)*Rn(:,4);
    BC.C(:,3)=(15/8)*Rn(:,2)-(35/4)*Rn(:,4)+(63/8)*Rn(:,6);
    
    %D1,D2,D3
    BC.D(:,1)=(3/2)*Rn(:,3)-(5/2)*Rn(:,5);
    BC.D(:,2)=(3/2)*Rn(:,3)-(5/2)*Rn(:,5);
    BC.D(:,3)=(-15/8)*Rn(:,3)+(35/4)*Rn(:,5)-(63/8)*Rn(:,7);
    
    %E1,E2,E3
    BC.E(:,1)=(15/8)*Rn(:,2)-(35/4)*Rn(:,4)+(63/8)*Rn(:,6);
    BC.E(:,2)=(-45/16)*Rn(:,2)+(285/16)*Rn(:,4)-(539/16)*Rn(:,6)+(315/16)*Rn(:,8);
    BC.E(:,3)=(-45/16)*Rn(:,2)+(285/16)*Rn(:,4)-(539/16)*Rn(:,6)+(315/16)*Rn(:,8);
    
    %F1,F2,F3
    BC.F(:,1)=(-15/8)*Rn(:,3)+(35/4)*Rn(:,5)-(63/8)*Rn(:,7);
    BC.F(:,2)=(45/16)*Rn(:,3)-(285/16)*Rn(:,5)+(539/16)*Rn(:,7)-(315/16)*Rn(:,9);
    BC.F(:,3)=(45/16)*Rn(:,3)-(285/16)*Rn(:,5)+(539/16)*Rn(:,7)-(315/16)*Rn(:,9);
    
    ksi=[BC.A BC.B BC.C BC.D BC.E BC.F];
 
    clear BC
    
elseif sp==7
    %A1,A2,A3,A4
    BC.A(:,1)=-Rn(:,2);
    BC.A(:,2)=(-9/4)*Rn(:,2)+(15/2)*Rn(:,4)-(25/4)*Rn(:,6);
    BC.A(:,3)=(-225/64)*Rn(:,2)+(525/16)*Rn(:,4)-(3395/32)*Rn(:,6)+...
        (2205/16)*Rn(:,8)-(3969/64)*Rn(:,10);
    BC.A(:,4)=(-1225/256)*Rn(:,2)+(11025/128)*Rn(:,4)-(147735/256)*Rn(:,6)+...
     (116655/64)*Rn(:,8)-(750519/256)*Rn(:,10)+(297297/128)*Rn(:,12)-...
     (184041/256)*Rn(:,14);
    
    %B1, B2, B3,B4
    BC.B(:,1)=3*Rn(:,3);
    BC.B(:,2)=(63/4)*Rn(:,3)-(105/2)*Rn(:,5)+(175/4)*Rn(:,7);
    BC.B(:,3)=(2475/64)*Rn(:,3)-(5775/16)*Rn(:,5)+(37345/32)*Rn(:,7)-...
        (24255/16)*Rn(:,9)+(43659/64)*Rn(:,11);
    BC.B(:,4)=15*((1225/256)*Rn(:,3)-(11025/128)*Rn(:,5)+(147735/256)*Rn(:,7)-...
     (116655/64)*Rn(:,9)+(750519/256)*Rn(:,11)-(297297/128)*Rn(:,13)+(184041/256)*Rn(:,15));
    
    %C1, C2,C3,C4
    BC.C(:,1)=(-3/2)*Rn(:,2)+(5/2)*Rn(:,4);
    BC.C(:,2)=(-3/2)*Rn(:,2)+(5/2)*Rn(:,4);
    BC.C(:,3)=(15/8)*Rn(:,2)-(35/4)*Rn(:,4)+(63/8)*Rn(:,6);
    BC.C(:,4)=(-35/16)*Rn(:,2)+(315/16)*Rn(:,4)-(693/16)*Rn(:,6)+(429/16)*Rn(:,8);
    
    %D1,D2,D3,D4
    BC.D(:,1)=(3/2)*Rn(:,3)-(5/2)*Rn(:,5);
    BC.D(:,2)=(3/2)*Rn(:,3)-(5/2)*Rn(:,5);
    BC.D(:,3)=(-15/8)*Rn(:,3)+(35/4)*Rn(:,5)-(63/8)*Rn(:,7);
    BC.D(:,4)=(35/16)*Rn(:,3)-(315/16)*Rn(:,5)+(693/16)*Rn(:,7)-(429/16)*Rn(:,9);
    
    %E1,E2,E3,E4
    BC.E(:,1)=(15/8)*Rn(:,2)-(35/4)*Rn(:,4)+(63/8)*Rn(:,6);
    BC.E(:,2)=(-45/16)*Rn(:,2)+(285/16)*Rn(:,4)-(539/16)*Rn(:,6)+(315/16)*Rn(:,8);
    BC.E(:,3)=(-45/16)*Rn(:,2)+(285/16)*Rn(:,4)-(539/16)*Rn(:,6)+(315/16)*Rn(:,8);
    BC.E(:,4)=(105/32)*Rn(:,2)-35*Rn(:,4)+(1827/16)*Rn(:,6)-(297/2)*Rn(:,8)+(2145/32)*Rn(:,10);
    
    %F1,F2,F3,F4
    BC.F(:,1)=(-15/8)*Rn(:,3)+(35/4)*Rn(:,5)-(63/8)*Rn(:,7);
    BC.F(:,2)=(45/16)*Rn(:,3)-(285/16)*Rn(:,5)+(539/16)*Rn(:,7)-(315/16)*Rn(:,9);
    BC.F(:,3)=(45/16)*Rn(:,3)-(285/16)*Rn(:,5)+(539/16)*Rn(:,7)-(315/16)*Rn(:,9);
    BC.F(:,4)=(-105/32)*Rn(:,3)+35*Rn(:,5)-(1827/16)*Rn(:,7)+(297/2)*Rn(:,9)-(2145/32)*Rn(:,11);
    
    %G1,G2,G3,G4
    BC.G(:,1)=(-35/16)*Rn(:,2)+(315/16)*Rn(:,4)-(693/16)*Rn(:,6)+(429/16)*Rn(:,8);
    BC.G(:,2)=(105/32)*Rn(:,2)-35*Rn(:,4)+(1827/16)*Rn(:,6)-(297/2)*Rn(:,8)+(2145/32)*Rn(:,10);
    BC.G(:,3)=(-525/128)*Rn(:,2)+(7175/128)*Rn(:,4)-(17325/64)*Rn(:,6)+(37395/64)*Rn(:,8)-...
    (73689/128)*Rn(:,10)+(27027/128)*Rn(:,12);
    BC.G(:,4)=(-525/128)*Rn(:,2)+(7175/128)*Rn(:,4)-(17325/64)*Rn(:,6)+(37395/64)*Rn(:,8)-...
    (73689/128)*Rn(:,10)+(27027/128)*Rn(:,12);
    
    %H1,H2,H3,H4
    BC.H(:,1)=(35/16)*Rn(:,3)-(315/16)*Rn(:,5)+(693/16)*Rn(:,7)-(429/16)*Rn(:,9);
    BC.H(:,2)=(-105/32)*Rn(:,3)+35*Rn(:,5)-(1827/16)*Rn(:,7)+(297/2)*Rn(:,9)-(2145/32)*Rn(:,11);
    BC.H(:,3)=(525/128)*Rn(:,3)-(7175/128)*Rn(:,5)+(17325/64)*Rn(:,7)-(37395/64)*Rn(:,9)+...
        (73689/128)*Rn(:,11)-(27027/128)*Rn(:,13);
    BC.H(:,4)=(525/128)*Rn(:,3)-(7175/128)*Rn(:,5)+(17325/64)*Rn(:,7)-(37395/64)*Rn(:,9)+...
        (73689/128)*Rn(:,11)-(27027/128)*Rn(:,13);
    
    ksi=[BC.A BC.B BC.C BC.D BC.E BC.F BC.G BC.H];
%    mysave('sp7.bc',foo);
    
    clear BC
else
    ksi=[];
end
